---
title: "Replication Guide"
subtitle: "Interest Group Persusian Experiment"
format: 
  html:
    number-sections: true
    embed-resources: true
    toc: true
    code-fold: true
---



# Replication and data 
All analyses were conducted in R, with package versions captured via `sessioninfo::session_info()`. The file `replication.qmd` reproduces every figure and table in the order they appear in the paper (main text and appendix). Random seeds are set for all stochastic procedures (`set.seed(201911)` for general analyses; `set.seed(3452892)` for forests). The code reads `1_data/df_all.rds` (analysis data) and `1_data/vars.xlsx` (variable labels) and writes outputs to `3_figures/` and `2_tables/`. A complete, runnable bundle (code, data, and outputs) is provided in the replication package; knitting `replication.qmd` recreates all results.


## R Session Info

```{r setup, echo=FALSE, message=FALSE,warning=FALSE}
rm(list=ls())

knitr::opts_chunk$set(echo=TRUE, warning=FALSE, message=FALSE, comment=NA)
options(qwraps2_markup = "markdown")

## Packages
if (!require(pacman)) install.packages("pacman")

pacman::p_load(
  ggpubr,     # plot arrange
  estimatr,
  gridExtra,
  modelsummary,
  #ggridges,    # plot
  readxl,      # read excel
  aod,         # hypothesis testing (1.3.1)
  car,         # linear hypothesis testing for causal tree (3.0-2)
  dplyr,       # Data manipulation (0.8.0.1)
  evtree,      # evolutionary learning of globally optimal trees (1.0-7)
  fastDummies,
  fBasics,     # Summary statistics (3042.89)
  ggplot2,
  grf,         # Generalized random forests (0.10.2)
  #haven,       # load sav
  kableExtra,  # Prettier RMarkdown (1.0.1)
  knitr,
  labelled,
  purrr,
  rpart,       # Classification and regression trees, or CART (4.1-13)
  reshape2,
  rpart.plot,  # Plotting trees (3.0.6)
  readr,       # Reading csv files (1.3.1)
  sjlabelled,
  stats,
  summarytools,
  texreg,
  tidyverse,
  tidyselect,
  tidyr,       # Database operations (0.8.3)
  treeClust,   # Predicting leaf position for causal trees (1.1-7)
  tibble,
  broom.helpers,
  naniar,   # create facror scores
  psych,     # create factor scores
  vtable,   # Summary Stats
  broom,
  ggtext,   # ggplot add on
  viridisLite,
  here,
  ggh4x,
  emmeans,
  patchwork
  )      

sessioninfo::session_info()


# Set seed for reproducibility
set.seed(201911) 

# Paths
# Paths
here("replication.qmd")

```


## Load data
```{r,message=FALSE,warning=FALSE}

# variable names and labels
var_list <- read_excel("1_data/vars.xlsx") 
# data
df <- readRDS("1_data/df_all.rds")


```

## Load functions
```{r ,message=FALSE,warning=FALSE}
# functions (to calculate factor scores)
save_fa_scores <- function(d){
  fa.mod <- fa(d,1)
  factor.scores(d,fa.mod)$scores
}

```

# Main Analysis

## Figure 1: Main Effects
```{r }

## helper: run lm_robust() and return a tidy row
make_summary <- function(df, y, coalition, strategy, policy, outcome) {
  df %>%
    group_by(
      coalition = .data[[coalition]],
      strategy  = .data[[strategy]]
    ) %>%
    do(tidy(lm_robust(reformulate("1", y), data = .))) %>%
    transmute(
      policy, outcome, coalition, strategy,
      estimate = estimate,
      conf.low, conf.high
    )
}

summary_long <- bind_rows(
  make_summary(df, "ecar_choice",  "ecar_coalition", "ecar_strategy",
               "E-car subsidies", "Choice"),
  make_summary(df, "co2_choice",   "co2_coalition",  "co2_strategy",
               " CO2 tax",         "Choice"),
  make_summary(df, "ecar_letter",  "ecar_coalition", "ecar_strategy",
               "E-car subsidies", "Letter"),
  make_summary(df, "co2_letter",   "co2_coalition",  "co2_strategy",
               " CO2 tax",         "Letter")
)


summary_long <- summary_long %>% 
  ## turn ecar_choice / co2_choice → “Choice”, the others → “Letter”
  mutate(
    outcome = case_when(
      str_detect(outcome, "choice")  ~ "Choice",
      str_detect(outcome, "letter")  ~ "Letter",
      TRUE                           ~ NA_character_
    ),
    ## make it an ordered factor (keeps strip order)
    outcome = factor(outcome, levels = c("Choice", "Letter"))
  )

fig_1 <- ggplot(summary_long,
              aes(coalition, estimate,
                  colour = strategy,      # keep colour legend
                  group  = strategy,
                  shape  = strategy,      # we still draw shapes …
                  linetype = strategy)) + # … and linetypes
  geom_point(aes(shape = strategy),
             position = position_dodge(0.5),
             size = 3, alpha = .4) +
  geom_line(aes(linetype = strategy),
            position = position_dodge(0.5)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                position = position_dodge(0.5), width = 0) +
  facet_grid(outcome ~ policy, scales = "free_y") +
  ## hide shape & linetype guides, keep colour
  guides(shape     = "none",
         linetype  = "none") +
  theme_bw(base_size = 14) +
  theme(axis.text.x     = element_text(angle = 45, hjust = 1),
        strip.text      = element_text(face = "bold"),
        legend.position = "bottom") +
  labs(x = "Support",
       y = "Estimate",
       colour = "Strategy")   # title for the *remaining* legend

fig_1
ggsave("3_figures/fig_1.pdf", width = 20, height = 20, units = "cm")

```


## Figure 2: Baseline support
```{r }


# co2
summary_df_co2 <-
  df %>%
  select(co2_coalition,co2_strategy,co2_choice) %>%
  filter(co2_strategy == "Text" & co2_coalition == "Control") %>%
  add_column(control = "co2") %>%
  select(control,co2_choice) %>%
  rename(outcome = co2_choice) 


# ecar
summary_df_ecar <-
  df %>%
  select(ecar_coalition,ecar_strategy,ecar_choice) %>%
  filter(ecar_strategy == "Text" & ecar_coalition == "Control") %>%
  add_column(control = "ecar") %>%
  select(control,ecar_choice) %>%
  rename(outcome = ecar_choice) 

control_df<-rbind(summary_df_co2,summary_df_ecar)


basline <- lh_robust(
  outcome ~ 0 + control,             # cell-means coding
  data   = control_df,
  linear_hypothesis = "controlco2 = controlecar")


# plot

basline_df<-tidy(basline) %>%
  filter(term %in% c("controlco2", "controlecar")) %>%     # keep only the two terms
  dplyr::mutate(term = dplyr::recode(term,                                # rename them
                       "controlco2"  = "CO2",
                       "controlecar" = "Ecar"))

  baseline <-
  ggplot(basline_df, aes(x=term, y=estimate)) +
  geom_point(position=position_dodge(width=0.5), alpha=0.4, size = 3) +
  geom_line(position=position_dodge(width=0.5)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), position=position_dodge(width=0.5), width = 0) +
  theme_bw() +
  theme(text = element_text(size=12))+
  ylab("Policy Support [1 = 'Yes', 0 = 'No']")+
  xlab("Control")+ 
  theme(legend.position="bottom") + 
  labs(color = "Type Support")+ 
  guides(group = FALSE, shape = FALSE,linetype = FALSE)+
  ggtitle("Baseline support")

  baseline
  
  ggsave("3_figures/fig_2.pdf", width = 15, height = 7.5, units = "cm")



```



## Figure 3: Main Effects by country

```{r }


make_summary <- function(df, y,
                         coalition, strategy,
                         policy, outcome) {

  df %>% 
    mutate(country = factor(cov_country,           # 0/1 → readable labels
                            levels   = c(0, 1),
                            labels   = c("Germany", "UK"))) %>% 
    group_by(country,
             coalition = .data[[coalition]],
             strategy  = .data[[strategy]]) %>% 
    do(tidy(lm_robust(reformulate("1", y), data = .))) %>% 
    transmute(country, policy, outcome,
              coalition, strategy,
              estimate = estimate,
              conf.low, conf.high)
}

summary_long <- bind_rows(
  make_summary(df, "ecar_choice", "ecar_coalition", "ecar_strategy",
               "E-car subsidies", "Choice"),
  make_summary(df, "co2_choice",  "co2_coalition",  "co2_strategy",
               "CO2 tax",         "Choice"),
  make_summary(df, "ecar_letter", "ecar_coalition", "ecar_strategy",
               "E-car subsidies", "Letter"),
  make_summary(df, "co2_letter",  "co2_coalition",  "co2_strategy",
               "CO2 tax",         "Letter")
)%>% mutate(outcome = case_when( str_detect(outcome, regex("choice", ignore_case = TRUE)) ~ "Choice", str_detect(outcome, regex("letter", ignore_case = TRUE)) ~ "Letter" ), outcome = factor(outcome, levels = c("Choice", "Letter")))

# ── build the plot exactly as before … ──────────────────────────────────────────
dodge <- position_dodge(0.5)

plt <- ggplot(summary_long,
              aes(coalition, estimate,
                  colour   = strategy,
                  group    = interaction(strategy, country),
                  shape    = country,
                  linetype = strategy)) +
  geom_point(position = dodge, size = 3, alpha = .4) +
  geom_line(position  = dodge) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                position = dodge, width = 0)+
  theme_bw()

# ── NEW: nested faceting ───────────────────────────────────────────────────────
plt <- plt +
  facet_nested(outcome  ~ policy+country,        # row strips: Country ▸ Outcome
               labeller = label_value,
               scales   = "free_y" ) +           # prints "country = X", "outcome = Y"
  theme(panel.spacing = unit(0, "line"),          # removes blank space between panels
        strip.text.y  = element_text(angle = 0),  # keep row-strip text horizontal
        axis.text.x   = element_text(angle = 45, hjust = 1),
        strip.text    = element_text(face = "bold"),
        legend.position = "bottom") +
  labs(x      = "Support",
       y      = "Estimate",
       colour = "Strategy",
       shape  = "Country") +
  guides(                         # 1 ─ hide unwanted guides
    shape    = "none",            # ← Country
    linetype = "none"             # ← strategy (dashes)
  ) +
  labs(colour = NULL)             # 2 ─ blank title (or use "")


print(plt)

ggsave("3_figures/fig_1_country.pdf", width = 15, height = 15, units = "cm")




```

## Table 4
Note that the regression underlying Table 4 is estimated on the same data we use for the causal forest; therefore, we run the causal forest first.

```{r }
# Causal Forests  =====
# Number of trees for causal forests

# Set seed for reproducibility
set.seed(3452892)

# Additional Functions
to_dummy <- function(x){
  ux <- unique(x)
  K <- length(ux)
  for(i in 1:length(ux)){
    x[x==ux[i]] <- i-1
  }
  return(as.numeric(x))
}

# prepare covariates
covariates<-df %>%
  dplyr::select(vars_select(names(df), starts_with('cov', ignore.case = TRUE)))    %>%
  #dplyr::select(-c(contains("_dk")))   %>%
  dplyr::select(-c(contains("position")))   %>%
  dplyr::select(-c(contains("industry")))   %>%
  dplyr::select(-c(cov_alignment_gp,
                   cov_alignment_indusry1,
                   cov_alignment_indusry2,
                   cov_alignment_union1,
                   cov_alignment_union2,
                   cov_alignment_climate_alliance,
                   cov_home_sqm,
                   cov_log_home_sqm))   


# Ensure covariates are a character vector
covariate_names <- colnames(covariates)
# prepare outcomes
df<-df  %>%
    dplyr::mutate(
    ecar_arg_factor = save_fa_scores(dplyr::select(.,update_support_ecar_emissions:update_support_ecar_effective)),
    ecar_decr_factor = save_fa_scores(dplyr::select(.,update_des_ecar,update_des_pers_ecar)),
    ecar_inj_factor = save_fa_scores(dplyr::select(.,update_inj_ecar,update_inj_pers_ecar)),
    )

# Define outcomes
out_list <- list(
  "ecar" = c(
    "ecar_arg_factor",
    "update_salience_ecar",
    "ecar_decr_factor",
    "ecar_inj_factor"
  )
)

# prepare binary treatments

df <- df %>%
  mutate(
    ecar_combined = paste0(ecar_coalition, "_", ecar_strategy),
    ecar_combined = case_when(
      ecar_coalition == "Control" & ecar_strategy == "Text" ~ "Control",
      TRUE ~ ecar_combined
    )
  )


df <- model.matrix(~ 0 + ecar_combined, df) %>%
  as.data.frame() %>%
  bind_cols(df, .)

treat_group_list <- unique(df$ecar_combined)
treat_group_list <- treat_group_list[treat_group_list != "Control"]


df_hte<-df %>% 
    dplyr::select(ID,ecar_coalition,ecar_strategy,ecar_choice,ecar_letter,
                  ecar_arg_factor,update_salience_ecar,ecar_decr_factor,
                  ecar_inj_factor) %>%
    dplyr::bind_cols(covariates) %>% 
    drop_na()

ecar.main<-lm_robust(ecar_choice ~ecar_strategy*ecar_coalition,data=df_hte) 
ecar.letter<-lm_robust(ecar_letter ~ecar_strategy*ecar_coalition,data=df_hte) 
ecar.ecar_arg_factor<-lm_robust(ecar_arg_factor ~ecar_strategy*ecar_coalition,data=df_hte) 
ecar.update_salience_ecar<-lm_robust(update_salience_ecar ~ecar_strategy*ecar_coalition,data=df_hte) 
ecar.ecar_decr_factor<-lm_robust(ecar_decr_factor ~ecar_strategy*ecar_coalition,data=df_hte) 
ecar.ecar_inj_factor<-lm_robust(ecar_inj_factor ~ecar_strategy*ecar_coalition,data=df_hte) 


hte.coef.map = list(
     "(Intercept)" = "Intercept", 
     "ecar_coalitionIG" = "IG", 
     "ecar_coalitionIG+Coalition" = "IG+Coalition", 
     "ecar_coalitionIG+Protest" = "IG+Protest", 
     "ecar_coalitionIG+Coalition+Protest" = "IG+Coalition+Protest"
     )


mod_names <- c("Choice",
               "Letter",
               "Argument",
               "Salience", 
               "Public Comp.",
               "Public Support")


table1_agg_hte<-texreg(list(ecar.main,
                          ecar.letter,
                          ecar.ecar_arg_factor,
                          ecar.update_salience_ecar,
                          ecar.ecar_decr_factor,
                          ecar.ecar_inj_factor
                          ),
       include.ci = FALSE,
       custom.model.names=mod_names,
       caption = "E-car Policy, based on fully saturated models. Results for interest group strategy are not displayed. Sample conisists of data with complete set of covariates also used for causal forests.",
       #omit.coef=c('(Intercept)'),
       stars = c(0.01, 0.05, 0.1),
       custom.coef.map = hte.coef.map,
       label = "table1_hte:coefficients"
       )

write(table1_agg_hte,file="2_tables/tab_4.tex")

html_code <- texreg::htmlreg(
  list(ecar.main, ecar.letter, ecar.ecar_arg_factor, ecar.update_salience_ecar,
       ecar.ecar_decr_factor, ecar.ecar_inj_factor),
  include.ci = FALSE,
  custom.model.names = mod_names,
  caption = "E-car Policy, based on fully saturated models. Results for interest group strategy are not displayed. Sample consists of data with a complete set of covariates also used for causal forests.",
  stars = c(0.01, 0.05, 0.1),
  custom.coef.map = hte.coef.map,
  doctype = FALSE
)

knitr::asis_output(html_code)  

```


# Appendix


## Summary Statistics



###  Summary statistics: Socioeconomics and Preferences
```{r , message=FALSE}


summ_stats <-df %>%
  dplyr::select(vars_select(names(df), starts_with('cov', ignore.case = TRUE)))    %>%
  dplyr::select(-c(contains("bin")))   %>%
  dplyr::select(-c(cov_valid)) %>%
  fBasics::basicStats() %>%
  t() %>%
  as.data.frame() %>%
  dplyr::select("NAs", "Mean", "Stdev", "Minimum", "Median", "Maximum") 

# Add in labels
summ_stats1 <- summ_stats %>% 
  mutate(Variable = factor(rownames(summ_stats), var_list$new_name, var_list$label)) %>% 
  relocate(Variable)  %>% 
  dplyr::filter( !grepl("Industry|Membership|Job|Alignment", Variable))

# Pretty-printing in HTML
summ_stats_table <- kable(summ_stats1, "html", digits = 2, booktabs = TRUE, row.names = FALSE)

kable_styling(summ_stats_table,
              bootstrap_options=c("striped", "hover", "condensed", "responsive"),
              full_width=FALSE)


tab_sum_stats1 <- kable(summ_stats1, format = "latex", digits = 2, caption = "Summary statistics: Socioeconomics and Preferences", booktabs = T, linesep = "", label = "SummStats", row.names = FALSE) %>%
  kable_styling(latex_options="scale_down")


write(tab_sum_stats1,file="2_tables/sum_stats1.tex")

```

###  Summary statistics: Employment and Occupation
```{r , message=FALSE}


summ_stats2 <- summ_stats %>% 
  mutate(Variable = factor(rownames(summ_stats), var_list$new_name, var_list$label)) %>% 
  relocate(Variable) %>% 
  dplyr::filter( grepl("Industry|Membership|Job", Variable))

tab_sum_stats2 <- kable(summ_stats2, "html", digits = 2, booktabs = TRUE, row.names = FALSE)
tab_sum_stats2


tab_sum_stats2 <- kable(summ_stats2, format = "latex", digits = 2, caption = "Summary statistics: Employment and Occupation", booktabs = T, linesep = "", label = "SummStats", row.names = FALSE) %>%
  kable_styling(latex_options="scale_down")


write(tab_sum_stats2,file="2_tables/sum_stats2.tex")


```

###  Summary statistics: Preferences

Note that the question on the alignment between respondent and position of organizations has many missing values. This is because many respondents answered "Dont know organisation".
```{r , message=FALSE}

summ_stats3 <- summ_stats %>% 
  mutate(Variable = factor(rownames(summ_stats), var_list$new_name, var_list$label)) %>% 
  relocate(Variable) %>% 
  dplyr::filter( grepl("Alignment", Variable))


tab_summ_stats3 <- kable(summ_stats3, "html", digits = 2, booktabs = TRUE, row.names = FALSE)
tab_summ_stats3


tab_sum_stats3 <- kable(summ_stats3, format = "latex", digits = 2, caption = "Summary statistics: Preferences", booktabs = T, linesep = "", label = "SummStats", row.names = FALSE) %>%
  kable_styling(latex_options="scale_down")


write(tab_sum_stats3,file="2_tables/sum_stats3.tex")
```


## Further results

### Main Effects

```{r, message=FALSE}
df <- readRDS("1_data/df_all.rds")

long.coef.map = list(
     "(Intercept)" = "Intercept", 
     "ecar_coalitionIG" = "Support IG", 
     "ecar_coalitionIG+Coalition" = "Support IG+Coalition", 
     "ecar_coalitionIG+Protest" = "Support IG+Protest", 
     "ecar_coalitionIG+Coalition+Protest" = "Support IG+Coalition+Protest", 
     "ecar_strategyFlyer" = "Support Flyer", 
     "ecar_strategyVideo" = "Support Video", 
     "ecar_strategyFlyer:ecar_coalitionIG" = "Support IG $\\times$ Support Flyer",
     "ecar_strategyVideo:ecar_coalitionIG" = "Support IG $\\times$ Support Video",
     "ecar_strategyFlyer:ecar_coalitionIG+Coalition" = "Support IG+Coalition  $\\times$ Support Flyer",
     "ecar_strategyVideo:ecar_coalitionIG+Coalition" = "Support IG+Coalition  $\\times$ Support Video",
     "ecar_strategyFlyer:ecar_coalitionIG+Coalition+Protest" = "Support IG+Coalition+Protest $\\times$ Support Flyer",
     "ecar_strategyVideo:ecar_coalitionIG+Coalition+Protest" = "Support IG+Coalition+Protest  $\\times$ Support Video",
     "ecar_strategyFlyer:ecar_coalitionIG+Protest" = "Support IG+Protest  $\\times$ Flyer",
     "ecar_strategyVideo:ecar_coalitionIG+Protest" = "Support IG+Protest  $\\times$ Video"
     )


ecar.main<-lm_robust(ecar_choice ~ecar_strategy*ecar_coalition,data=df) 
ecar.letter<-lm_robust(ecar_letter ~ecar_strategy*ecar_coalition,data=df) 



table1_long<-texreg(list(ecar.main,ecar.letter),
       include.ci = FALSE,
       custom.model.names=c("Support","Letter"),
       caption = "E-car Policy, fully saturated models",
       #omit.coef=c('(Intercept)'),
       stars = c(0.01, 0.05, 0.1),
       custom.coef.map = long.coef.map,
       label = "table1_long:coefficients"
       )

write(table1_long,file="2_tables/tab_D1_long.tex")


html_code <- texreg::htmlreg(
  list(ecar.main,ecar.letter),
  custom.model.names=c("Support","Letter"),
  caption = "E-car Policy, fully saturated models",
  custom.coef.map = long.coef.map,
  include.ci = FALSE,
  doctype = FALSE
)

knitr::asis_output(html_code)  

```



```{r , message=FALSE}
# create index
df<-df  %>%
    dplyr::mutate(
    co2_arg_fa = save_fa_scores(dplyr::select(.,update_CO2taxation_compensation:update_CO2taxation_research)),
    co2_decr_fa = save_fa_scores(dplyr::select(.,update_des_co2,update_des_pers_co2)),
    co2_inj_fa = save_fa_scores(dplyr::select(.,update_inj_co2,update_inj_pers_co2))
    )


long.coef.map = list(
     "(Intercept)" = "Intercept", 
     "co2_coalitionIG" = "Support IG", 
     "co2_coalitionIG+Coalition" = "Support IG+Coalition", 
     "co2_coalitionIG+Coalition+Protest" = "Support IG+Coalition+Protest", 
     "co2_strategyFlyer" = "Support Flyer", 
     "co2_strategyVideo" = "Support Video", 
     "co2_strategyFlyer:co2_coalitionIG" = "Support IG $\\times$ Support Flyer",
     "co2_strategyVideo:co2_coalitionIG" = "Support IG $\\times$ Support Video",
     "co2_strategyFlyer:co2_coalitionIG+Coalition" = "Support IG+Coalition  $\\times$ Support Flyer",
     "co2_strategyVideo:co2_coalitionIG+Coalition" = "Support IG+Coalition  $\\times$ Support Video",
     "co2_strategyFlyer:co2_coalitionIG+Protest" = "Support IG+Protest $\\times$  Support Flyer",
     "co2_strategyVideo:co2_coalitionIG+Protest" = "Support IG+Protest  $\\times$  Support Video",
     "co2_strategyFlyer:co2_coalitionIG+Coalition+Protest" = "Support IG+Coalition+Protest $\\times$ Support Flyer",
     "co2_strategyVideo:co2_coalitionIG+Coalition+Protest" = "Support IG+Coalition+Protest  $\\times$ Support Video"
     )



co2.choice<-lm_robust(co2_choice ~co2_strategy*co2_coalition,data=df) 
co2.letter<-lm_robust(co2_letter ~co2_strategy*co2_coalition,data=df) 


table2_long<-texreg(list(co2.choice,co2.letter),
       include.ci = FALSE,
       custom.model.names=c("Support","Letter"),
       caption = "CO2 Policy, fully saturated models",
       stars = c(0.01, 0.05, 0.1),
       custom.coef.map = long.coef.map,
       label = "table2_long:coefficients"
       )

write(table2_long,file="2_tables/tab_D2_long.tex.tex")

html_code <- texreg::htmlreg(
 list(co2.choice,co2.letter),
       include.ci = FALSE,
       custom.model.names=c("Support","Letter"),
       caption = "CO2 Policy, fully saturated models",
       stars = c(0.01, 0.05, 0.1),
       custom.coef.map = long.coef.map,
        doctype = FALSE
)

knitr::asis_output(html_code)  


```


### Hypothesis H1a: Strategy

```{r , message=FALSE}

## Fit the factorial model once -----------------------------
fit <- lm_robust(ecar_choice ~ ecar_coalition * ecar_strategy,
                 data = df, se_type = "stata")


# Specify the four coalition cues we want
camp_levels <- c("IG", "IG+Coalition", "IG+Protest", "IG+Coalition+Protest")

# 2. Get campaign-only marginal means of strategy
emm_strat_camp <- emmeans(
  fit,
  ~ ecar_strategy,
  at       = list(ecar_coalition = camp_levels),  # <-- excludes Control
  weights  = "equal"                              # or "proportional"
)

# 3. Pairwise or planned contrasts
h1a_camp <- contrast(
  emm_strat_camp,
  method = list("Video – Text" = c(-1, 0, 1),
                "Flyer – Text" = c(-1, 1, 0)),
  adjust = "none")

summary(h1a_camp, infer = c(TRUE, TRUE))


if (knitr::is_latex_output()) {
  # Save LaTeX table (no visible output in the doc)
  kbl(h1a_camp, format = "latex", booktabs = TRUE,
      caption = "Contrasts testing H1a, Ecar: Video/Flyer vs Text") |>
    kable_styling(latex_options = c("hold_position")) |>
    save_kable("2_tables/tab_D3.tex")
}

if (knitr::is_html_output()) {
  # Show HTML table in HTML output
  kbl(h1a_camp, format = "html", booktabs = TRUE,
      caption = "Contrasts testing H1a, Ecar: Video/Flyer vs Text")
}
```


```{r , message=FALSE}


## Fit the factorial model once -----------------------------
fit <- lm_robust(co2_choice ~ co2_coalition * co2_strategy,
                 data = df, se_type = "stata")


# Specify the four coalition cues we want
camp_levels <- c("IG", "IG+Coalition", "IG+Protest", "IG+Coalition+Protest")

# 2. Get campaign-only marginal means of strategy
emm_strat_camp <- emmeans(
  fit,
  ~ co2_strategy,
  at       = list(co2_coalition = camp_levels),  # <-- excludes Control
  weights  = "equal"                              # or "proportional"
)

# 3. Pairwise or planned contrasts
h1a_camp <- contrast(
  emm_strat_camp,
  method = list("Video – Text" = c(-1, 0, 1),
                "Flyer – Text" = c(-1, 1, 0)),
  adjust = "none")

summary(h1a_camp, infer = c(TRUE, TRUE))


if (knitr::is_latex_output()) {
  # Save LaTeX table (no visible output in PDF)
  kbl(h1a_camp, format = "latex", booktabs = TRUE,
      caption = "Contrasts testing H1a, CO2: Video/Flyer vs Text") |>
    kable_styling(latex_options = "hold_position") |>
    save_kable("2_tables/tab_D4.tex")
}

if (knitr::is_html_output()) {
  # Show HTML table in HTML output
  kbl(h1a_camp, format = "html", booktabs = TRUE,
      caption = "Contrasts testing H1a, CO2: Video/Flyer vs Text")
}
```



### Linear Hypothesis (H1b, H2a, H2b)

```{r , message=FALSE}

# Flyer vs. control
H1b1<-lh_robust(ecar_choice ~ ecar_coalition*ecar_strategy, data = df, linear_hypothesis = "ecar_strategyFlyer=0") %>% tidy() %>% filter(str_detect(term, '=')) %>%  add_column(variable = "H1b: Strategy Flyer vs. control") 
# Video vs. control
H1b2<-lh_robust(ecar_choice ~ ecar_coalition*ecar_strategy, data = df, linear_hypothesis = "ecar_strategyVideo=0") %>% tidy() %>% filter(str_detect(term, '=')) %>%  add_column(variable = "H1b: Strategy Video vs. control") 
# IG coalition vs. control
H2a1<-lh_robust(ecar_choice ~ ecar_coalition*ecar_strategy, data = df, linear_hypothesis = "ecar_coalitionIG+Coalition=0") %>% tidy() %>% filter(str_detect(term, '=')) %>%  add_column(variable = "H2a: IG+coalition vs. control")  
# IG protest vs. control
H2a2<-lh_robust(ecar_choice ~ ecar_coalition*ecar_strategy, data = df, linear_hypothesis = "ecar_coalitionIG+Protest=0") %>% tidy() %>% filter(str_detect(term, '=')) %>%  add_column(variable = "H2a: IG+protest vs. control")  
# IG coalition vs. IG
H2b1<-lh_robust(ecar_choice ~ ecar_coalition*ecar_strategy, data = df, linear_hypothesis = "ecar_coalitionIG = ecar_coalitionIG+Coalition") %>% tidy() %>% filter(str_detect(term, '=')) %>%  add_column(variable = "H2b: IG+coalition vs. IG")  
# IG coalition vs. IG
H2b2<-lh_robust(ecar_choice ~ ecar_coalition*ecar_strategy, data = df, linear_hypothesis = "ecar_coalitionIG = ecar_coalitionIG+Protest") %>% tidy() %>% filter(str_detect(term, '=')) %>%  add_column(variable = "H2b: IG+protest vs. IG")  

# bind
hyp<-rbind(H1b1,H1b2,H2a1,H2a2,H2b1,H2b2) %>% mutate(outcome = ifelse(outcome == "co2_choice", "Support", NA))

#create texreg
outcome <- unique(as.character(hyp$outcome))
tr <- lapply(outcome, function(x) {
  d <- hyp[hyp$variable == x, ]
  t <- createTexreg(coef.names = as.character(hyp$variable),
                    coef = hyp$estimate,
                    se = hyp$std.error,
                    pvalues = hyp$p.value,
                    model.name = x)
  return(t)
})

lh_ecar<-texreg(tr,caption = "Ecar Policy Support, Linear hypothesis based on fully saturated model")

write(lh_ecar,file="2_tables/tab_D5.tex")

html_code <- texreg::htmlreg(tr,
       caption = "Ecar Policy Support, fully saturated models"
)

knitr::asis_output(html_code)  


```



```{r , message=FALSE}

# Flyer vs. control
H1b1<-lh_robust(co2_choice ~ co2_coalition*co2_strategy, data = df, linear_hypothesis = "co2_strategyFlyer=0") %>% tidy() %>% filter(str_detect(term, '=')) %>%  add_column(variable = "H1b: Strategy Flyer vs. control") 
# Video vs. control
H1b2<-lh_robust(co2_choice ~ co2_coalition*co2_strategy, data = df, linear_hypothesis = "co2_strategyVideo=0") %>% tidy() %>% filter(str_detect(term, '=')) %>%  add_column(variable = "H1b: Strategy Video vs. control") 
# IG coalition vs. control
H2a1<-lh_robust(co2_choice ~ co2_coalition*co2_strategy, data = df, linear_hypothesis = "co2_coalitionIG+Coalition=0") %>% tidy() %>% filter(str_detect(term, '=')) %>%  add_column(variable = "H2a: IG+coalition vs. control")  
# IG protest vs. control
H2a2<-lh_robust(co2_choice ~ co2_coalition*co2_strategy, data = df, linear_hypothesis = "co2_coalitionIG+Protest=0") %>% tidy() %>% filter(str_detect(term, '=')) %>%  add_column(variable = "H2a: IG+protest vs. control")  
# IG coalition vs. IG
H2b1<-lh_robust(co2_choice ~ co2_coalition*co2_strategy, data = df, linear_hypothesis = "co2_coalitionIG = co2_coalitionIG+Coalition") %>% tidy() %>% filter(str_detect(term, '=')) %>%  add_column(variable = "H2b: IG+coalition vs. IG")  
# IG coalition vs. IG
H2b2<-lh_robust(co2_choice ~ co2_coalition*co2_strategy, data = df, linear_hypothesis = "co2_coalitionIG = co2_coalitionIG+Protest") %>% tidy() %>% filter(str_detect(term, '=')) %>%  add_column(variable = "H2b: IG+protest vs. IG")  

# bind
hyp<-rbind(H1b1,H1b2,H2a1,H2a2,H2b1,H2b2) %>% mutate(outcome = ifelse(outcome == "co2_choice", "Support", NA))

#create texreg
outcome <- unique(as.character(hyp$outcome))
tr <- lapply(outcome, function(x) {
  d <- hyp[hyp$variable == x, ]
  t <- createTexreg(coef.names = as.character(hyp$variable),
                    coef = hyp$estimate,
                    se = hyp$std.error,
                    pvalues = hyp$p.value,
                    model.name = x)
  return(t)
})

lh_co2<-texreg(tr,caption = "CO2 Policy Support, Linear hypothesis based on fully saturated model")

write(lh_co2,file="2_tables/tab_D6.tex")


html_code <- texreg::htmlreg(tr,
       caption = "CO2 Policy, fully saturated models"
)

knitr::asis_output(html_code)  

```


### Linear Hypothesis, Baseline (H3)

```{r , message=FALSE}

# co2
summary_df_co2 <-
  df %>%
  select(co2_coalition,co2_strategy,co2_choice) %>%
  filter(co2_strategy == "Text" & co2_coalition == "Control") %>%
  add_column(control = "co2") %>%
  select(control,co2_choice) %>%
  rename(outcome = co2_choice) 


# ecar
summary_df_ecar <-
  df %>%
  select(ecar_coalition,ecar_strategy,ecar_choice) %>%
  filter(ecar_strategy == "Text" & ecar_coalition == "Control") %>%
  add_column(control = "ecar") %>%
  select(control,ecar_choice) %>%
  rename(outcome = ecar_choice) 

control_df<-rbind(summary_df_co2,summary_df_ecar)


basline <- lh_robust(
  outcome ~ 0 + control,             # cell-means coding
  data   = control_df,
  linear_hypothesis = "controlco2 = controlecar")


# table
baseline.coef.map = list(
     "controlco2" = "Control CO2", 
     "controlecar" = "Control Ecar",
     "controlco2 = controlecar"= "Control CO2 = Control Ecar"
     )

basline_tab<-texreg(basline,
       include.ci = FALSE,
       caption = "Baseline support, E-car and CO2 Policies",
       stars = c(0.01, 0.05, 0.1),
       custom.coef.map = baseline.coef.map,
       label = "tab_baseline"
       )


write(basline_tab,file="2_tables/tab_D7.tex")

html_code <- texreg::htmlreg(basline,
       caption = "Baseline support, E-car and CO2 Policies",
       stars = c(0.01, 0.05, 0.1),
       custom.coef.map = baseline.coef.map)

knitr::asis_output(html_code)  


```


### Further results: Mechanisms


## cates treatment - outcome
```{r }

# Causal Forests  =====
# Number of trees for causal forests

# Set seed for reproducibility
set.seed(3452892)


# Define outcomes
out_list <- list(
  "ecar" = c(
    "ecar_choice",
    "ecar_letter"
)
)



# prepare outcomes
df<-df  %>%
    dplyr::mutate(
    ecar_arg_factor = save_fa_scores(dplyr::select(.,update_support_ecar_emissions:update_support_ecar_effective)),
    ecar_decr_factor = save_fa_scores(dplyr::select(.,update_des_ecar,update_des_pers_ecar)),
    ecar_inj_factor = save_fa_scores(dplyr::select(.,update_inj_ecar,update_inj_pers_ecar)),
    )

# also select mediator for cate

ecar.main<-lm_robust(ecar_choice ~ecar_strategy*ecar_coalition,data=df_hte) 
ecar.letter<-lm_robust(ecar_letter ~ecar_strategy*ecar_coalition,data=df_hte) 
ecar.ecar_arg_factor<-lm_robust(ecar_arg_factor ~ecar_strategy*ecar_coalition,data=df_hte) 
ecar.update_salience_ecar<-lm_robust(update_salience_ecar ~ecar_strategy*ecar_coalition,data=df_hte) 
ecar.ecar_decr_factor<-lm_robust(ecar_decr_factor ~ecar_strategy*ecar_coalition,data=df_hte) 
ecar.ecar_inj_factor<-lm_robust(ecar_inj_factor ~ecar_strategy*ecar_coalition,data=df_hte) 

# prepare covariates
covariates<-df %>%
  dplyr::select(c(
    vars_select(names(df), starts_with('cov', ignore.case = TRUE)),
    "ecar_arg_factor",
    "update_salience_ecar",
    "ecar_decr_factor",
    "ecar_inj_factor"
  )) %>%
  #dplyr::select(-c(contains("_dk")))   %>%
  dplyr::select(-c(contains("position")))   %>%
  dplyr::select(-c(contains("industry")))   %>%
  dplyr::select(-c(cov_alignment_gp,
                   cov_alignment_indusry1,
                   cov_alignment_indusry2,
                   cov_alignment_union1,
                   cov_alignment_union2,
                   cov_alignment_climate_alliance,
                   cov_home_sqm,
                   cov_log_home_sqm))   


# Ensure covariates are a character vector
covariate_names <- colnames(covariates)

# prepare binary treatments
df <- df %>%
  mutate(
    ecar_combined = paste0(ecar_coalition, "_", ecar_strategy),
    ecar_combined = case_when(
      ecar_coalition == "Control" & ecar_strategy == "Text" ~ "Control",
      TRUE ~ ecar_combined
    )
  )


df <- model.matrix(~ 0 + ecar_combined, df) %>%
  as.data.frame() %>%
  bind_cols(df, .)

treat_group_list <- unique(df$ecar_combined)
treat_group_list <- treat_group_list[treat_group_list != "Control"]


# run forests
out_cat <- 1
out <- 1
N_trees = 2000
train_fraction <- 1  

# Initialize a dataframe to store all CATEs
cates_all <- data.frame()

# Iterate over outcome categories and outcomes
for (out_cat in 1:length(out_list)) {
  outcomes <- out_list[[out_cat]]
  
  for (out in 1:length(outcomes)) {
    cat("Running causal forests for", outcomes[out], "\n")
    
    # Iterate over treatment groups
    for (i in 1:length(treat_group_list)) {
      # Filter dataset to include control and current treatment
      df_cf <- df %>%
        filter(ecar_combined %in% c("Control", treat_group_list[i])) %>%
        select(all_of(c("ecar_combined", outcomes[out], covariate_names))) %>%
        drop_na()
      
      # Convert covariates to numeric dummies
      df_cf <- df_cf %>%
        mutate_at(all_of(covariate_names), to_dummy) %>%
        mutate(
          ecar_combined = ifelse(ecar_combined == "Control", 0, 1),
          !!outcomes[out] := as.numeric(.[[outcomes[out]]])
        )
      
      # Fit the causal forest
      c.forest <- causal_forest(
        X = as.matrix(df_cf[, covariate_names]),
        Y = df_cf[[outcomes[out]]],
        W = df_cf$ecar_combined,
        num.trees = N_trees
      )
      
      # Predict CATEs
      oob_pred <- predict(c.forest, estimate.variance = TRUE)
      oob_tauhat_cf <- oob_pred$predictions
      oob_tauhat_cf_se <- sqrt(oob_pred$variance.estimates)
      
      # Create dataframe for the current treatment and outcome
      df_oob_cate <- data.frame(
        cate = oob_tauhat_cf,
        cate_se = oob_tauhat_cf_se,
        treat = treat_group_list[i],
        outcome = outcomes[out]
      )
      
      # Combine into the main CATEs dataframe
      cates_all <- bind_rows(cates_all, df_oob_cate)
    }
  }
}

 
# Update coalition categories
cates_all <- cates_all %>%
  mutate(
    strategy = case_when(
      grepl("Flyer", treat) ~ "Flyer",
      grepl("Text", treat) ~ "Text",
      grepl("Video", treat) ~ "Video",
      TRUE ~ "Other"
    ),
    outcome = case_when(
      outcome == "ecar_choice" ~ "Choice",
      outcome == "ecar_letter" ~ "Letter",
      TRUE ~ outcome  # Default: Keep the original name if no match
        ),
    coalition = case_when(
      grepl("^Control", treat) ~ "Control",  # Treatments starting with "Control"
      grepl("^IG_Flyer|^IG_Text|^IG_Video", treat) ~ "IG",  # IG alone without coalition
      grepl("IG\\+Coalition\\+Protest", treat) ~ "IG+Coalition+Protest",  # Explicit match for "IG+Coalition+Protest"
      grepl("IG\\+Coalition", treat) ~ "IG+Coalition",  # Explicit match for "IG+Coalition"
      grepl("IG\\+Protest", treat) ~ "IG+Protest",  # Explicit match for "IG+Protest"
      TRUE ~ "Other"  # Catch-all for unexpected cases
    )
  )


ggplot(cates_all, aes(x = cate, fill = strategy, color = strategy)) +
  geom_histogram(alpha = 0.4, bins = 30, position = "identity") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  labs(
    x = "Estimated CATE",
    y = "Frequency",
    title = "Distribution of Estimated CATEs of Treatment on Outcomes"
  ) +
  facet_grid(coalition ~ outcome, scales = "free") +  # Facet by coalition and outcome
  scale_fill_manual(values = c(
    "Flyer" = "#1f77b4",
    "Text" = "#ff7f0e",
    "Video" = "#2ca02c"
  )) +
  scale_color_manual(values = c(
    "Flyer" = "#1f77b4",
    "Text" = "#ff7f0e",
    "Video" = "#2ca02c"
  )) +
  theme_bw() +
  theme(
    strip.text = element_text(size = 10, face = "bold"),
    strip.background = element_rect(fill = "lightgray"),
    legend.position = "bottom",
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 12, face = "bold")
  )


ggsave("3_figures/fig_D1.pdf", width = 23, height = 23, units = "cm")


```


## Figure 4 causal forest
```{r,warning = FALSE }

# Causal Forests  =====
# Number of trees for causal forests

# Set seed for reproducibility
set.seed(3452892)

# Additional Functions
to_dummy <- function(x){
  ux <- unique(x)
  K <- length(ux)
  for(i in 1:length(ux)){
    x[x==ux[i]] <- i-1
  }
  return(as.numeric(x))
}

# prepare covariates
covariates<-df %>%
  dplyr::select(vars_select(names(df), starts_with('cov', ignore.case = TRUE)))    %>%
  #dplyr::select(-c(contains("_dk")))   %>%
  dplyr::select(-c(contains("position")))   %>%
  dplyr::select(-c(contains("industry")))   %>%
  dplyr::select(-c(cov_alignment_gp,
                   cov_alignment_indusry1,
                   cov_alignment_indusry2,
                   cov_alignment_union1,
                   cov_alignment_union2,
                   cov_alignment_climate_alliance,
                   cov_home_sqm,
                   cov_log_home_sqm))   


# Ensure covariates are a character vector
covariate_names <- colnames(covariates)
# prepare outcomes
df<-df  %>%
    dplyr::mutate(
    ecar_arg_factor = save_fa_scores(dplyr::select(.,update_support_ecar_emissions:update_support_ecar_effective)),
    ecar_decr_factor = save_fa_scores(dplyr::select(.,update_des_ecar,update_des_pers_ecar)),
    ecar_inj_factor = save_fa_scores(dplyr::select(.,update_inj_ecar,update_inj_pers_ecar)),
    )

# Define outcomes
out_list <- list(
  "ecar" = c(
    "ecar_arg_factor",
    "update_salience_ecar",
    "ecar_decr_factor",
    "ecar_inj_factor"
  )
)

# prepare binary treatments

df <- df %>%
  mutate(
    ecar_combined = paste0(ecar_coalition, "_", ecar_strategy),
    ecar_combined = case_when(
      ecar_coalition == "Control" & ecar_strategy == "Text" ~ "Control",
      TRUE ~ ecar_combined
    )
  )


df <- model.matrix(~ 0 + ecar_combined, df) %>%
  as.data.frame() %>%
  bind_cols(df, .)

treat_group_list <- unique(df$ecar_combined)
treat_group_list <- treat_group_list[treat_group_list != "Control"]


# run forests
out_cat <- 1
out <- 1
N_trees = 2000
train_fraction <- 1  

# Initialize a dataframe to store all CATEs
cates_all <- data.frame()

# Iterate over outcome categories and outcomes
for (out_cat in 1:length(out_list)) {
  outcomes <- out_list[[out_cat]]
  
  for (out in 1:length(outcomes)) {
    cat("Running causal forests for", outcomes[out], "\n")
    
    # Iterate over treatment groups
    for (i in 1:length(treat_group_list)) {
      # Filter dataset to include control and current treatment
      df_cf <- df %>%
        filter(ecar_combined %in% c("Control", treat_group_list[i])) %>%
        select(all_of(c("ecar_combined", outcomes[out], covariate_names))) %>%
        drop_na()
      
      # Convert covariates to numeric dummies
      df_cf <- df_cf %>%
        mutate_at(all_of(covariate_names), to_dummy) %>%
        mutate(
          ecar_combined = ifelse(ecar_combined == "Control", 0, 1),
          !!outcomes[out] := as.numeric(.[[outcomes[out]]])
        )
      
      # Fit the causal forest
      c.forest <- causal_forest(
        X = as.matrix(df_cf[, covariate_names]),
        Y = df_cf[[outcomes[out]]],
        W = df_cf$ecar_combined,
        num.trees = N_trees
      )
      
      # Predict CATEs
      oob_pred <- predict(c.forest, estimate.variance = TRUE)
      oob_tauhat_cf <- oob_pred$predictions
      oob_tauhat_cf_se <- sqrt(oob_pred$variance.estimates)
      
      # Create dataframe for the current treatment and outcome
      df_oob_cate <- data.frame(
        cate = oob_tauhat_cf,
        cate_se = oob_tauhat_cf_se,
        treat = treat_group_list[i],
        outcome = outcomes[out]
      )
      
      # Combine into the main CATEs dataframe
      cates_all <- bind_rows(cates_all, df_oob_cate)
    }
  }
}


# Update coalition categories
cates_all <- cates_all %>%
  mutate(
    strategy = case_when(
      grepl("Flyer", treat) ~ "Flyer",
      grepl("Text", treat) ~ "Text",
      grepl("Video", treat) ~ "Video",
      TRUE ~ "Other"
    ),
    outcome = case_when(
      outcome == "ecar_arg_factor" ~ "Argument",
      outcome == "update_salience_ecar" ~ "Salience",
      outcome == "ecar_decr_factor" ~ "Public Compliance",
      outcome == "ecar_inj_factor" ~ "Public Support",
      TRUE ~ outcome  # Default: Keep the original name if no match
        ),
    coalition = case_when(
      grepl("^Control", treat) ~ "Control",  # Treatments starting with "Control"
      grepl("^IG_Flyer|^IG_Text|^IG_Video", treat) ~ "IG",  # IG alone without coalition
      grepl("IG\\+Coalition\\+Protest", treat) ~ "IG+Coalition+Protest",  # Explicit match for "IG+Coalition+Protest"
      grepl("IG\\+Coalition", treat) ~ "IG+Coalition",  # Explicit match for "IG+Coalition"
      grepl("IG\\+Protest", treat) ~ "IG+Protest",  # Explicit match for "IG+Protest"
      TRUE ~ "Other"  # Catch-all for unexpected cases
    )
  )


ggplot(cates_all, aes(x = cate, fill = strategy, color = strategy)) +
  geom_histogram(alpha = 0.4, bins = 30, position = "identity") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  labs(
    x = "Estimated CATE",
    y = "Frequency",
    title = "Distribution of Estimated CATEs for the effect the treatment on the mediotors"
  ) +
  facet_grid(coalition ~ outcome, scales = "free") +  # Facet by coalition and outcome
  scale_fill_manual(values = c(
    "Flyer" = "#1f77b4",
    "Text" = "#ff7f0e",
    "Video" = "#2ca02c"
  )) +
  scale_color_manual(values = c(
    "Flyer" = "#1f77b4",
    "Text" = "#ff7f0e",
    "Video" = "#2ca02c"
  )) +
  theme_bw() +
  theme(
    strip.text = element_text(size = 10, face = "bold"),
    strip.background = element_rect(fill = "lightgray"),
    legend.position = "bottom",
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 12, face = "bold")
  )


ggsave("3_figures/fig_D2.pdf", width = 23, height = 23, units = "cm")

```


### Mechanisms: Disaggregated

```{r, message = FALSE,warning=FALSE}
# Causal Forests  =====
# Number of trees for causal forests

# Set seed for reproducibility
set.seed(3452892)

# Define outcomes
out_list <- list(
  "ecar" = c(
    "update_support_ecar_emissions",
    "update_support_ecar_noise",
    "update_support_ecar_cooperation",
    "update_support_ecar_costs",
    "update_support_ecar_effective"
)
)


# prepare covariates
covariates<-df %>%
  dplyr::select(c(
    vars_select(names(df), starts_with('cov', ignore.case = TRUE)),
    "ecar_arg_factor",
    "update_salience_ecar",
    "ecar_decr_factor",
    "ecar_inj_factor"
  )) %>%
  #dplyr::select(-c(contains("_dk")))   %>%
  dplyr::select(-c(contains("position")))   %>%
  dplyr::select(-c(contains("industry")))   %>%
  dplyr::select(-c(cov_alignment_gp,
                   cov_alignment_indusry1,
                   cov_alignment_indusry2,
                   cov_alignment_union1,
                   cov_alignment_union2,
                   cov_alignment_climate_alliance,
                   cov_home_sqm,
                   cov_log_home_sqm))   


# Ensure covariates are a character vector
covariate_names <- colnames(covariates)

# prepare binary treatments
df <- df %>%
  mutate(
    ecar_combined = paste0(ecar_coalition, "_", ecar_strategy),
    ecar_combined = case_when(
      ecar_coalition == "Control" & ecar_strategy == "Text" ~ "Control",
      TRUE ~ ecar_combined
    )
  )


df <- model.matrix(~ 0 + ecar_combined, df) %>%
  as.data.frame() %>%
  bind_cols(df, .)

treat_group_list <- unique(df$ecar_combined)
treat_group_list <- treat_group_list[treat_group_list != "Control"]


# run forests
out_cat <- 1
out <- 1
N_trees = 2000
train_fraction <- 1  

# Initialize a dataframe to store all CATEs
cates_all <- data.frame()

# Iterate over outcome categories and outcomes
for (out_cat in 1:length(out_list)) {
  outcomes <- out_list[[out_cat]]
  
  for (out in 1:length(outcomes)) {
    cat("Running causal forests for", outcomes[out], "\n")
    
    # Iterate over treatment groups
    for (i in 1:length(treat_group_list)) {
      # Filter dataset to include control and current treatment
      df_cf <- df %>%
        filter(ecar_combined %in% c("Control", treat_group_list[i])) %>%
        select(all_of(c("ecar_combined", outcomes[out], covariate_names))) %>%
        drop_na()
      
      # Convert covariates to numeric dummies
      df_cf <- df_cf %>%
        mutate_at(all_of(covariate_names), to_dummy) %>%
        mutate(
          ecar_combined = ifelse(ecar_combined == "Control", 0, 1),
          !!outcomes[out] := as.numeric(.[[outcomes[out]]])
        )
      
      # Fit the causal forest
      c.forest <- causal_forest(
        X = as.matrix(df_cf[, covariate_names]),
        Y = df_cf[[outcomes[out]]],
        W = df_cf$ecar_combined,
        num.trees = N_trees
      )
      
      # Predict CATEs
      oob_pred <- predict(c.forest, estimate.variance = TRUE)
      oob_tauhat_cf <- oob_pred$predictions
      oob_tauhat_cf_se <- sqrt(oob_pred$variance.estimates)
      
      # Create dataframe for the current treatment and outcome
      df_oob_cate <- data.frame(
        cate = oob_tauhat_cf,
        cate_se = oob_tauhat_cf_se,
        treat = treat_group_list[i],
        outcome = outcomes[out]
      )
      
      # Combine into the main CATEs dataframe
      cates_all <- bind_rows(cates_all, df_oob_cate)
    }
  }
}

 
# Update coalition categories
cates_all <- cates_all %>%
  mutate(
    strategy = case_when(
      grepl("Flyer", treat) ~ "Flyer",
      grepl("Text", treat) ~ "Text",
      grepl("Video", treat) ~ "Video",
      TRUE ~ "Other"
    ),
    outcome = case_when(
      outcome == "update_support_ecar_emissions" ~ "Emissions",
      outcome == "update_support_ecar_noise" ~ "Noise",
      outcome == "update_support_ecar_cooperation" ~ "Cooperation",
      outcome == "update_support_ecar_costs" ~ "Costs",
      outcome == "update_support_ecar_effective" ~ "Effective",
      TRUE ~ outcome  # Default: Keep the original name if no match
        ),
    coalition = case_when(
      grepl("^Control", treat) ~ "Control",  # Treatments starting with "Control"
      grepl("^IG_Flyer|^IG_Text|^IG_Video", treat) ~ "IG",  # IG alone without coalition
      grepl("IG\\+Coalition\\+Protest", treat) ~ "IG+Coalition+Protest",  # Explicit match for "IG+Coalition+Protest"
      grepl("IG\\+Coalition", treat) ~ "IG+Coalition",  # Explicit match for "IG+Coalition"
      grepl("IG\\+Protest", treat) ~ "IG+Protest",  # Explicit match for "IG+Protest"
      TRUE ~ "Other"  # Catch-all for unexpected cases
    )
  )



ggplot(cates_all, aes(x = cate, fill = strategy, color = strategy)) +
  geom_histogram(alpha = 0.4, bins = 30, position = "identity") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  labs(
    x = "Estimated CATE",
    y = "Frequency",
    title = "Distribution of Estimated CATEs of Treatment on Argument Mediators"
  ) +
  facet_grid(coalition ~ outcome, scales = "free") +  # Facet by coalition and outcome
  scale_fill_manual(values = c(
    "Flyer" = "#1f77b4",
    "Text" = "#ff7f0e",
    "Video" = "#2ca02c"
  )) +
  scale_color_manual(values = c(
    "Flyer" = "#1f77b4",
    "Text" = "#ff7f0e",
    "Video" = "#2ca02c"
  )) +
  theme_bw() +
  theme(
    strip.text = element_text(size = 10, face = "bold"),
    strip.background = element_rect(fill = "lightgray"),
    legend.position = "bottom",
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 12, face = "bold")
  )


ggsave("3_figures/fig_D3.pdf", width = 23, height = 23, units = "cm")



```


```{r}

# Causal Forests  =====
# Number of trees for causal forests

# Set seed for reproducibility
set.seed(3452892)


# Define outcomes
out_list <- list(
  "ecar" = c(
    "update_inj_ecar",
    "update_inj_pers_ecar"
)
)


# prepare covariates
covariates<-df %>%
  dplyr::select(c(
    vars_select(names(df), starts_with('cov', ignore.case = TRUE)),
    "ecar_arg_factor",
    "update_salience_ecar",
    "ecar_decr_factor",
    "ecar_inj_factor"
  )) %>%
  #dplyr::select(-c(contains("_dk")))   %>%
  dplyr::select(-c(contains("position")))   %>%
  dplyr::select(-c(contains("industry")))   %>%
  dplyr::select(-c(cov_alignment_gp,
                   cov_alignment_indusry1,
                   cov_alignment_indusry2,
                   cov_alignment_union1,
                   cov_alignment_union2,
                   cov_alignment_climate_alliance,
                   cov_home_sqm,
                   cov_log_home_sqm))   


# Ensure covariates are a character vector
covariate_names <- colnames(covariates)

# prepare binary treatments
df <- df %>%
  mutate(
    ecar_combined = paste0(ecar_coalition, "_", ecar_strategy),
    ecar_combined = case_when(
      ecar_coalition == "Control" & ecar_strategy == "Text" ~ "Control",
      TRUE ~ ecar_combined
    )
  )


df <- model.matrix(~ 0 + ecar_combined, df) %>%
  as.data.frame() %>%
  bind_cols(df, .)

treat_group_list <- unique(df$ecar_combined)
treat_group_list <- treat_group_list[treat_group_list != "Control"]


# run forests
out_cat <- 1
out <- 1
N_trees = 2000
train_fraction <- 1  

# Initialize a dataframe to store all CATEs
cates_all <- data.frame()

# Iterate over outcome categories and outcomes
for (out_cat in 1:length(out_list)) {
  outcomes <- out_list[[out_cat]]
  
  for (out in 1:length(outcomes)) {
    cat("Running causal forests for", outcomes[out], "\n")
    
    # Iterate over treatment groups
    for (i in 1:length(treat_group_list)) {
      # Filter dataset to include control and current treatment
      df_cf <- df %>%
        filter(ecar_combined %in% c("Control", treat_group_list[i])) %>%
        select(all_of(c("ecar_combined", outcomes[out], covariate_names))) %>%
        drop_na()
      
      # Convert covariates to numeric dummies
      df_cf <- df_cf %>%
        mutate_at(all_of(covariate_names), to_dummy) %>%
        mutate(
          ecar_combined = ifelse(ecar_combined == "Control", 0, 1),
          !!outcomes[out] := as.numeric(.[[outcomes[out]]])
        )
      
      # Fit the causal forest
      c.forest <- causal_forest(
        X = as.matrix(df_cf[, covariate_names]),
        Y = df_cf[[outcomes[out]]],
        W = df_cf$ecar_combined,
        num.trees = N_trees
      )
      
      # Predict CATEs
      oob_pred <- predict(c.forest, estimate.variance = TRUE)
      oob_tauhat_cf <- oob_pred$predictions
      oob_tauhat_cf_se <- sqrt(oob_pred$variance.estimates)
      
      # Create dataframe for the current treatment and outcome
      df_oob_cate <- data.frame(
        cate = oob_tauhat_cf,
        cate_se = oob_tauhat_cf_se,
        treat = treat_group_list[i],
        outcome = outcomes[out]
      )
      
      # Combine into the main CATEs dataframe
      cates_all <- bind_rows(cates_all, df_oob_cate)
    }
  }
}

 
# Update coalition categories
cates_all <- cates_all %>%
  mutate(
    strategy = case_when(
      grepl("Flyer", treat) ~ "Flyer",
      grepl("Text", treat) ~ "Text",
      grepl("Video", treat) ~ "Video",
      TRUE ~ "Other"
    ),
    outcome = case_when(
      outcome == "update_inj_ecar" ~ "Public supprt, General Population",
      outcome == "update_inj_pers_ecar" ~ "Public supprt, Personal Environment",
      TRUE ~ outcome  # Default: Keep the original name if no match
        ),
    coalition = case_when(
      grepl("^Control", treat) ~ "Control",  # Treatments starting with "Control"
      grepl("^IG_Flyer|^IG_Text|^IG_Video", treat) ~ "IG",  # IG alone without coalition
      grepl("IG\\+Coalition\\+Protest", treat) ~ "IG+Coalition+Protest",  # Explicit match for "IG+Coalition+Protest"
      grepl("IG\\+Coalition", treat) ~ "IG+Coalition",  # Explicit match for "IG+Coalition"
      grepl("IG\\+Protest", treat) ~ "IG+Protest",  # Explicit match for "IG+Protest"
      TRUE ~ "Other"  # Catch-all for unexpected cases
    )
  )



ggplot(cates_all, aes(x = cate, fill = strategy, color = strategy)) +
  geom_histogram(alpha = 0.4, bins = 30, position = "identity") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  labs(
    x = "Estimated CATE",
    y = "Frequency",
    title = "Distribution of Estimated CATEs of Treatment on Mediators"
  ) +
  facet_grid(coalition ~ outcome, scales = "free") +  # Facet by coalition and outcome
  scale_fill_manual(values = c(
    "Flyer" = "#1f77b4",
    "Text" = "#ff7f0e",
    "Video" = "#2ca02c"
  )) +
  scale_color_manual(values = c(
    "Flyer" = "#1f77b4",
    "Text" = "#ff7f0e",
    "Video" = "#2ca02c"
  )) +
  theme_bw() +
  theme(
    strip.text = element_text(size = 10, face = "bold"),
    strip.background = element_rect(fill = "lightgray"),
    legend.position = "bottom",
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 12, face = "bold")
  )


ggsave("3_figures/fig_D4.pdf", width = 23, height = 23, units = "cm")


```

## Heterogeneous Treatment Effects

```{r, message = FALSE}
# get cov

covariates<-df %>%
  dplyr::select(vars_select(names(df), starts_with('cov', ignore.case = TRUE)))    %>%
  dplyr::select(-c(contains("bin")))   %>%
  dplyr::select(-c(contains("ind")))   %>%
  dplyr::select(-c(contains("dk")))   %>%
  #dplyr::select(-c(contains("alignment")))   %>%
  dplyr::select(-c(cov_valid,cov_age,cov_home_sqm))  


# get cov names
covariate_names <-colnames(covariates) 



# define binary treatments
df<-df %>% 
  dplyr::mutate(
      ecar_strategy_bin = dplyr::case_when(
                           ecar_strategy_num == 0 ~ 0,
                           ecar_strategy_num >=1 ~ 1),
      ecar_coalition_bin = dplyr::case_when(
                           ecar_coalition_num < 1 ~ 0,
                           ecar_coalition_num >=1 ~ 1))



# get data without NA's
df_hte<-df %>% 
    dplyr::select(ID,ecar_strategy_bin,ecar_coalition_bin,ecar_choice) %>% 
    dplyr::bind_cols(covariates) %>% 
    drop_na()

x<-df_hte %>% 
    summarise_all(typeof) %>% 
     gather

# Training fraction
train_fraction <- 1  # for causal forests: currently all data used


# A function to prep het effects data for a given treatment, removing others

# Labels for treatments
treatment_levels_diff <- c("ecar_strategy_bin","ecar_coalition_bin")
treatment_labels <- c("E-car strategy", "E-car Coalition")

num_tiles <- 4  # ntiles = CATE is above / below the median

# Note this is set up so that we can swap in the different outcome easily

het_df <- function(treatment_name = "ecar_strategy_bin", 
                   outcome_name = "ecar_choice",
                   data = df) {
  data$W = data[treatment_name][[1]]
  data$Y = data[outcome_name][[1]]
  data %>% 
    dplyr::select(Y, W, ID, all_of(covariate_names)) %>%
    drop_na() %>%
    mutate_if(is.factor, as.numeric)  %>%
    # mutate(federal.state =  as.numeric(factor(federal.state))) %>%
    # Trick to render W binary
    dplyr::mutate(
      Y = ifelse(W<0, -Y, Y),
      W = W^2)  
}




f_cf <- function(df, n_trees = 10000) {
  
  # rule of thumb: num.trees = number of individuals
  
  df_train <- sample_frac(df, replace=FALSE, size=train_fraction)
  df_test  <- anti_join(df, df_train, by = "ID") #need to check on larger 
  X <- as.matrix(df_train[, covariate_names])
  
  cf <-
    causal_forest(
    X = X,
    Y = df_train$Y,
    W = df_train$W,
    num.trees=n_trees) 

  #### Predict point estimates and standard errors (training set, out-of-bag)
  oob_pred      <- predict(cf, estimate.variance=TRUE)
  oob_tauhat_cf <- oob_pred$predictions
  oob_tauhat_cf_se <- sqrt(oob_pred$variance.estimates)
  
  #### Predict point estimates and standard errors (test set)
  # test_pred <- predict(cf, newdata=as.matrix(df_test[covariate_names]), estimate.variance=TRUE)
  # tauhat_cf_test <- test_pred$predictions
  # tauhat_cf_test_se <- sqrt(test_pred$variance.estimates)
  
  var_imp        <- c(variable_importance(cf)) 
  names(var_imp) <- covariate_names
  var_imp <- var_imp %>% sort(decreasing=TRUE)
  
  df_train$cate  <- oob_tauhat_cf
  df_train$ntile <- factor(ntile(oob_tauhat_cf, n=num_tiles))

# Standard model estimates by quantile
estimated_sample_ate <- 
  lm_robust(Y ~ ntile + ntile:W, data=df_train) %>% 
  tidy() %>% 
  dplyr::filter(stringr::str_detect(term, ":W"))

# AIPW estimates by quantile
estimated_aipw_ate <- 
  lapply(
  seq(num_tiles), function(w) 
    average_treatment_effect(cf, subset = df_train$ntile == w)
  ) %>% bind_rows

combined_estimates <- 
  bind_rows(
    estimated_sample_ate %>% mutate(type = "lm_robust") %>% dplyr::select(-outcome, -df, -statistic, - p.value),
    estimated_aipw_ate %>% rename(std.error = std.err) %>% 
      mutate(
        type  = "aipw",
        term = estimated_sample_ate$term,
        conf.low = estimate - 1.96*std.error,
        conf.high = estimate + 1.96*std.error)
  )
# Outputs
list(cf = cf,
     df_train = df_train, 
     df_test = df_test, 
     X = X,
     oob_tauhat_cf = oob_tauhat_cf, 
     var_imp = var_imp, 
     ntile_estimates = combined_estimates)
}


# Function to generate fitted values:
fitted_vals <- function(var_of_interest, model = test){
  
  df_train <- model$df_train
  cf <- model$cf
  
      is_continuous <- (length(unique(df_train[var_of_interest][[1]])) > 5) # crude rule for determining continuity
    if(is_continuous) {
      x_grid <- quantile(df_train[var_of_interest][[1]], probs = seq(0, 1, length.out = 5))
    } else {
      x_grid <- sort(unique(df_train[var_of_interest][[1]]))
    }
    
  df_grid <-  setNames(data.frame(x_grid), var_of_interest)
  
  other_covariates <- covariate_names[!covariate_names %in% var_of_interest]
  df_median <- df_train %>% dplyr::select(all_of(other_covariates)) %>% summarise_all(median) 
  df_eval <- crossing(df_median, df_grid)
  
  pred <- predict(cf, newdata=df_eval[,covariate_names], estimate.variance=TRUE)
df_eval$tauhat <- pred$predictions
df_eval$se <- sqrt(pred$variance.estimates)

# Change to factor so the plotted values are evenly spaced (e.g. logicals)
df_eval %>% arrange(var_of_interest) %>%
  mutate(var_of_interest = as.factor(as.numeric(df_eval[var_of_interest][[1]])))
}


cf_ecar <- 
  lapply(treatment_levels_diff, function(x) het_df(x, outcome_name = "ecar_choice") %>% f_cf(n_trees=10000))
names(cf_ecar) <- treatment_levels_diff

#Most common predictors

lapply(cf_ecar, function(j) names(j$var_imp[1:4])) %>% bind_rows(.id = "treatment") %>% kable(caption = "Strongest predictions", booktabs = TRUE)


what_matters_ecar <- 
  list(All = cf_ecar) %>%
  lapply(function(res) 
    lapply(res, function(j) j$var_imp %>% t %>% data.frame) %>%
      bind_rows %>% mutate(treatment = names(cf_ecar)) %>%
      gather(covariate, "value", - treatment)) %>% bind_rows(.id = "group")
  

what_matters_plot_ecar_top <- 
  what_matters_ecar  %>% 
  group_by(treatment) %>%
  arrange(value, .by_group = TRUE) %>%
  top_n(12) %>% 
  group_by(covariate) %>% 
  filter(n()>1)  %>% 
  mutate(treatment = case_when(treatment == 'ecar_strategy_bin' ~ 'strategy treatment',
                                 treatment == 'ecar_coalition_bin' ~ 'coalition treatment',
                                 TRUE ~ NA))    %>%
    dplyr::left_join(var_list, by=c('covariate'='new_name')) %>%
    dplyr::mutate(label = fct_rev(fct_reorder2(label,value,treatment))) %>%
    ggplot(aes(reorder(label, + value) , x=value, color=family)) +
    geom_point()+
    geom_segment(
    aes(y = label, yend = label, x = 0, xend = value, color=family)
    ) +    
    scale_y_discrete() +
    theme_bw() +
    facet_wrap(~treatment) +
    ylab(" ")+
    xlab("Variable Importance")+ 
    scale_x_continuous(limits = c(0,0.1 ))+
    theme(legend.position="bottom")+
    ggtitle("E-car")


##############################################
# CO2
##############################################

# define binary treatments
df<-df %>% 
  dplyr::mutate(
      co2_strategy_bin = dplyr::case_when(
                           co2_strategy_num == 0 ~ 0,
                           co2_strategy_num >=1 ~ 1),
      co2_coalition_bin = dplyr::case_when(
                           co2_coalition_num < 1 ~ 0,
                           co2_coalition_num >=1 ~ 1))


# get data without NA's
df_hte<-df %>% 
    dplyr::select(ID,co2_strategy_bin,co2_coalition_bin,co2_choice) %>% 
    dplyr::bind_cols(covariates) %>% 
    drop_na()

x<-df_hte %>% 
    summarise_all(typeof) %>% 
     gather


# Labels for treatments
treatment_levels_diff <- c("co2_strategy_bin","co2_coalition_bin")
treatment_labels <- c("CO2 strategy", "CO2 Coalition")

num_tiles <- 4  # ntiles = CATE is above / below the median

# Note this is set up so that we can swap in the different outcome easily

het_df <- function(treatment_name = "co2_strategy_bin", 
                   outcome_name = "co2_choice",
                   data = df) {
  data$W = data[treatment_name][[1]]
  data$Y = data[outcome_name][[1]]
  data %>% 
    dplyr::select(Y, W, ID, all_of(covariate_names)) %>%
    drop_na() %>%
    mutate_if(is.factor, as.numeric)  %>%
    # mutate(federal.state =  as.numeric(factor(federal.state))) %>%
    # Trick to render W binary
    dplyr::mutate(
      Y = ifelse(W<0, -Y, Y),
      W = W^2)  
}


# Run forests

cf_co2 <- 
  lapply(treatment_levels_diff, function(x) het_df(x, outcome_name = "co2_choice") %>% f_cf(n_trees=10000))
names(cf_co2) <- treatment_levels_diff

# Most common predictors

lapply(cf_co2, function(j) names(j$var_imp[1:6])) %>% bind_rows(.id = "treatment") %>% kable(caption = "Strongest predictions", booktabs = TRUE)


what_matters_co2<- 
  list(All = cf_co2) %>%
  lapply(function(res) 
    lapply(res, function(j) j$var_imp %>% t %>% data.frame) %>%
      bind_rows %>% mutate(treatment = names(cf_co2)) %>%
      gather(covariate, "value", - treatment)) %>% bind_rows(.id = "group")

what_matters_plot_co2_top <- 
  what_matters_co2  %>% 
  group_by(treatment) %>%
  arrange(value, .by_group = TRUE) %>%
  top_n(12) %>% 
  group_by(covariate) %>% 
  filter(n()>1)  %>% 
  dplyr::mutate(treatment = case_when(treatment == 'co2_strategy_bin' ~ 'strategy treatment',
                                 treatment == 'co2_coalition_bin' ~ 'coalition treatment',
                                 TRUE ~ NA))    %>%
    dplyr::left_join(var_list, by=c('covariate'='new_name')) %>%
    dplyr::mutate(label = fct_rev(fct_reorder2(label,value,treatment))) %>%
    ggplot(aes(reorder(label, + value) , x=value, color=family)) +
    geom_point()+
    geom_segment(
    aes(y = label, yend = label, x = 0, xend = value, color=family)
    ) +    
    scale_y_discrete() +
    theme_bw() +
    facet_wrap(~treatment) +
    ylab(" ")+
    xlab("Variable Importance")+ 
    scale_x_continuous(limits = c(0,0.1 ))+
    theme(legend.position="bottom")+
    ggtitle("CO2")



what_matters_plot_ecar_top + what_matters_plot_co2_top +
  plot_layout(ncol = 1, guides = "auto") &
  theme(legend.position = "bottom")


ggsave("3_figures/fig_D5.pdf", width = 10, height = 6, units = "in")

```


```{r, message = FALSE}
library("interplot")

# age
M1 <- lm(ecar_choice~cov_age_prec*ecar_strategy_bin*ecar_coalition_bin, data = df)
ecar_age  <- list(
  stratgey = interplot(m = M1, var1 = "ecar_strategy_bin", var2 = "cov_age_prec", plot = FALSE),
  colation = interplot(m = M1, var1 = "ecar_coalition_bin", var2 = "cov_age_prec", plot = FALSE)
  ) %>% dplyr::bind_rows(.id = "Treatment")  %>% mutate(variable="Age") %>% dplyr::rename(cov=cov_age_prec) %>%
ggplot(aes(x=cov, y=coef, colour=Treatment)) +
  geom_line() +
  geom_ribbon(aes(ymin=lb, ymax=ub, fill=Treatment), alpha=0.1) + 
  geom_hline(yintercept=0, linetype="longdash", lwd=0.35, size=0.75, colour = "#B55555") + 
  facet_wrap(~variable) +
  ylab("treatment effect") + 
  xlab("covariate")  + 
  theme_bw()+ theme(legend.position = "none")    

# sqm
M2 <- lm(ecar_choice~cov_log_home_sqm*ecar_strategy_bin*ecar_coalition_bin, data = df)
  ecar_sqm <- list(
  stratgey = interplot(m = M2, var1 = "ecar_strategy_bin", var2 = "cov_log_home_sqm", plot = FALSE),
  colation = interplot(m = M2, var1 = "ecar_coalition_bin", var2 = "cov_log_home_sqm", plot = FALSE)
  ) %>% bind_rows(.id = "Treatment") %>% mutate(variable="Home sqm (log)") %>% dplyr::rename(cov=cov_log_home_sqm)%>%
ggplot(aes(x=cov, y=coef, colour=Treatment)) +
  geom_line() +
  geom_ribbon(aes(ymin=lb, ymax=ub, fill=Treatment), alpha=0.1) + 
  geom_hline(yintercept=0, linetype="longdash", lwd=0.35, size=0.75, colour = "#B55555") + 
  facet_wrap(~variable) +
  ylab("treatment effect") + 
  xlab("covariate")  + 
  theme_bw()+ theme(legend.position = "none")  
  
# age
M1 <- lm(co2_choice~cov_age_prec*co2_strategy_bin*co2_coalition_bin, data = df)
age_co2 <- list(
  stratgey = interplot(m = M1, var1 = "co2_strategy_bin", var2 = "cov_age_prec", plot = FALSE),
  colation = interplot(m = M1, var1 = "co2_coalition_bin", var2 = "cov_age_prec", plot = FALSE)
  ) %>% dplyr::bind_rows(.id = "Treatment")  %>% mutate(variable="Age") %>% dplyr::rename(cov=cov_age_prec) %>%
ggplot(aes(x=cov, y=coef, colour=Treatment)) +
  geom_line() +
  geom_ribbon(aes(ymin=lb, ymax=ub, fill=Treatment), alpha=0.1) + 
  geom_hline(yintercept=0, linetype="longdash", lwd=0.35, size=0.75, colour = "#B55555") + 
  facet_wrap(~variable) +
  ylab("treatment effect") + 
  xlab("covariate")  + 
  theme_bw()

# sqm
M2 <- lm(co2_choice~cov_log_home_sqm*co2_strategy_bin*co2_coalition_bin, data = df)
sqm_co2 <- list(
  stratgey = interplot(m = M2, var1 = "co2_strategy_bin", var2 = "cov_log_home_sqm", plot = FALSE),
  colation = interplot(m = M2, var1 = "co2_coalition_bin", var2 = "cov_log_home_sqm", plot = FALSE)
  ) %>% bind_rows(.id = "Treatment") %>% mutate(variable="Home sqm (log)") %>% dplyr::rename(cov=cov_log_home_sqm)%>%
ggplot(aes(x=cov, y=coef, colour=Treatment)) +
  geom_line() +
  geom_ribbon(aes(ymin=lb, ymax=ub, fill=Treatment), alpha=0.1) + 
  geom_hline(yintercept=0, linetype="longdash", lwd=0.35, size=0.75, colour = "#B55555") + 
  facet_wrap(~variable) +
  ylab("treatment effect") + 
  xlab("covariate")  + 
  theme_bw()



co2_cov_plot<-ggarrange(age_co2,sqm_co2,
                         ncol=2,nrow=1,
                         common.legend = TRUE, legend="bottom")


ann1 <- ggplot() + 
  geom_text(aes(x=0, y=0, label = "bold(CO2)"), 
            parse = TRUE, size = 6, vjust = -.5) +
  theme_void()

ann2 <- ggplot() + 
  geom_text(aes(x=0, y=0, label = "bold(E-car)"), 
            parse = TRUE, size = 6, vjust = -0.5) +
  theme_void()
ecar_cov_plot<-ggarrange(ann2,ecar_sqm,ecar_age,
                         ann1,sqm_co2,age_co2,
                         ncol=3,nrow=2,
                         common.legend = TRUE, legend="bottom",
                         widths = c(0.2, 0.3, 0.3, 0.3)
)
ecar_cov_plot
ggsave("3_figures/fig_D6.pdf", width = 9, height = 5, units = "in")

```
